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ABSTRACT 

Heat  injur) >  is  a  problem  for  the  Armed  Forces,  especially  during  deployment  to  localities  with  very  hot  and 
humid  climates.  Early  warning  of  a  rising  core  body  temperature  (Tc)  can  help  prevent  heat  injuries.  To  this 
end,  we  developed  an  algorithm  that,  given  a  series  of  past  T(  measurements  obtained  using  an  digestible 
temperature  pill,  accumulates  evidence  of  a  rising  Tc  over  time  and  provides  ahead-of-time  warning  of  an 
impending,  dangerously  elevated  Tc.  Using  data  from  a  cohort  of  six  Soldiers  involved  in  field  exercises 
whose  Tc  exceeded  38.5  °C,  we  assessed  the  performance  of  the  warning  algorithm.  The  algorithm  predicted 
rises  in  Tc  with  a  clinically  useful  lead  time  (>  18  min)  and  reasonable  sensitivity  and  specificity  (>  87%). 
However,  because  digestible  temperature  pills  are  impractical  for  monitoring  a  large  number  of  Warfighters 
during  prolonged  operations,  we  developed  a  mathematical  model  that  uses  non-invasive  measurements  of 
physiological  variables,  such  as  activity  (A  c),  heart  rate  (Hr),  and  skin  temperature,  as  well  as 
environmental  information  [ambient  temperature  (TA)  and  relative  humidity  din)],  to  provide  individualized 
real-time  Tc  estimates.  Using  the  same  cohort  of  Soldiers,  we  evaluated  two  variants  of  the  individualized 
model,  one  that  used  all  the  measurements  (original  model)  and  another  that  used  A  c,  Hr,  and  month- 
average  Ta  and  /in  values  (reduced  model).  The  original  and  reduced  models  yielded  Tc  estimates  with 
average  errors  of  0.31X3  and  0.35  °C,  respectively,  which  are  within  the  physiological  intra-subject 
variability  of  Tc.  In  addition,  on  average,  the  estimation  time  delays  ranged  from  3  min  for  the  original 
model  to  5  min  for  the  reduced  model.  We  conclude  that  the  individualized  Tc  estimation  model  can  be  used 
to  replace  digestible  temperature  pills  and  enable  the  development  of  a  field-deployable  early  warning 
system  of  an  impending  rise  in  core  temperature. 

1.0  INTRODUCTION 

Heat-related  illnesses  are  a  significant  problem  for  the  United  States  (U.S.)  Armed  Forces,  especially  during 
deployments  to  localities  with  hot  and  humid  climates.  In  fact,  between  2009  and  2013  there  were  12,907 
heat  injuries  across  the  Services,  including  1,757  cases  of  heat  stroke  [1,  2].  The  primary  sign  of  an 
impending  heat  injury  is  a  rise  in  core  body  temperature  (Tc).  However,  during  strenuous  military  operations 
in  hot  and  humid  conditions,  Warfighters  focused  on  their  mission  could  miss  warning  signs  of  an  increasing 
Tc  and  an  impending  heat  illness  [3].  New  sensor  technologies,  which  afford  the  ability  to  measure  Tc  via 
ingestible  temperature  pills  [4],  and  mathematical  predictive  models  [5,  6]  could  be  coupled  to  develop  a 
hardware/software  warning  system  of  an  impending  rise  in  Tc  and  generate  alerts  to  potentially  prevent  heat 
injuries. 
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Although  ingestible  temperature  pills  provide  accurate  Tc  measurements,  they  are  invasive  and  cannot  be 
used  for  continuously  monitoring  large  number  of  Warfighters  during  prolonged  operations.  This  has  led  to 
the  investigation  of  more  practical  alternatives,  such  as  the  use  of  non-invasive  performance-  and  fitness¬ 
tracking  devices  that  monitor  a  wide  range  of  physiological  variables,  e.g.,  activity  (Ac),  heart  rate  (HR),  and 
skin  temperature  (Ts),  which  can  then  be  used  to  indirectly  infer  Tc.  This  could  be  accomplished  by  a 
physiology-based  mathematical  model  that  combines  the  information  from  such  non-invasive  measurements 
using  phenomenological  and  first-principles  concepts  to  provide  individualized  Tc  estimates  in  real  time. 

In  this  study,  we  report  on  the  development  of  algorithms  to  address  the  two  requirements  for  developing  a 
reliable  early  warning  system  for  heat  injury  prevention:  1)  an  algorithm  that,  given  a  time  series  of  Tc 
measurements  or  Tc  estimates,  provides  ahead-of-time  alerts  about  an  impending  rise  in  Tc  and  2)  an 
individualized  model  that  uses  non-invasive  measurements  of  Ac,  HR,  and  Ts  as  well  as  two  environmental 
variables  [ambient  temperature  ( TA )  and  relative  humidity  (Rh)]  to  provide  real-time  Tc  estimates. 

2.0  Tc  PREDICTION  AND  ALERT  ALGORITHMS 

Here,  we  detail  the  development  of  an  algorithm  that  uses  a  time  series  of  recent-past  Tc  measurements  to 
provide  reliable,  ahead-of-time  alerts  of  an  impending  rise  in  Tc. 

2.1  Methods:  Tc  Prediction  and  Alert  Algorithms 

2.1.1  Study  Data 

To  demonstrate  the  performance  of  the  proposed  algorithm,  we  used  data  from  a  field  study  involving  six 
U.S.  Army  Soldiers  [average  age:  23.1  year  (standard  deviation  [SD]  4.1);  average  height:  178  cm  (SD  7); 
average  weight:  81.3  kg  (SD  11.1)]  who  performed  regularly  scheduled  infantry  training.  The  training 
included  a  six-mile  foot  march  while  wearing  a  backpack  and  carrying  equipment  weighing  on  average  14.0 
kg  (SD  1.4)  and  exercises,  such  as  digging  ditches,  setting-up  concertina  wire,  performing  marksmanship 
drills,  running,  rolling,  and  jumping  as  part  of  the  approach  to  a  target.  During  the  training,  which  lasted  for  8 
to  14  h,  Soldiers  wore  the  advanced  combat  uniform  with  a  thermal  insulation  of  1.08  clo  and  an  evaporative 
potential  of  0.41  im/clo.  During  the  training  period,  the  recorded  average  TA  was  29.3°C  (SD  2.1),  the 
average  RH  was  74%  (SD  11),  and  the  average  wind  speed  was  2.5  m/s  (SD  0.8).  All  training  was  at  the 
direction  of  the  military  unit,  i.e.,  the  research  team  did  not  interfere  with  or  ask  for  any  alteration  to  training 
events.  The  Institutional  Review  Board  of  the  U.S.  Army  Research  Institute  of  Environmental  Medicine 
(Natick,  MA)  approved  the  study. 

Each  Soldier  was  instrumented  with  radio-thermometer  pills  (MiniMitter,  Inc.,  Bend,  OR)  that  measured  and 
transmitted  core  temperature  data  to  the  Hidalgo  Equivital  EQ-02  (Hidalgo,  Ltd.,  Cambridge,  UK) 
physiological  status  monitoring  (PSM)  system  [7],  Data  were  retrieved  from  the  PSM  system  at  the  end  of 
the  exercise  and  subsequently  analyzed.  Soldiers  ingested  the  pills  at  least  12  h  before  data  collection.  The 
thermometer  pills  had  the  following  technical  characteristics:  length:  21.9  mm;  diameter:  8.5  mm;  weight: 
1.75  grams;  sampling  period:  15  s;  temperature  range  32°C  -  42°C,  with  an  accuracy  of  ±0.1°C;  transmission 
method:  near-field  magnetic  link.  Additionally,  the  PSM  system  measured  the  following  three  physiological 
variables:  1 )  Ac  via  triaxial  accelerometers,  2)  HR  via  two  electrocardiogram  electrodes,  and  3 )  Ts  via  an 
infrared  sensor. 

2.1.2  Auto-regressive  (AR)  Model 

Given  temperature  measurements  TCf  .  at  discrete  time  t  sampled  every  S  min,  where  i  =  0,  1,. . .,  m- 1,  the  AR 
model  of  order  m  predicts  Tc  at  time  point  t+\,  TCt+v  through  a  linear  combination  of  the  antecedent  Tc 
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Figure  1 :  The  sequential  probability  ratio  test  (SPRT)  algorithm  that  provides  alerts  of  an  impending  rise  in 
core  body  temperature  (Tc).  The  auto-regressive  (AR)  model  uses  a  history  of  recent  Tc  measurements  to 
make  a  P-min-ahead  predictions.  It  also  provides  the  corresponding  prediction  intervals  (Pis),  i.e.,  the 
uncertainty  of  the  predictions.  The  SPRT  algorithm  uses  the  AR  model  outputs  to  decide  whether  Tc  would 
rise  beyond  a  predefined  temperature  threshold  and  provides  ahead-of-time  alerts  if  it  does. 

samples  and  AR  model  coefficients  b  (see  Figure  1).  To  make  predictions  M  time-steps  ahead  (prediction 
window  P  =  M  x  S  min  ahead),  we  iteratively  used  the  AR  model  equation  M  times,  substituting  the 
unobserved  signals  at  t  >  t+\  by  their  corresponding  predicted  values.  In  this  study,  we  set  the  sampling 
period  5  =  5  min  to  preserve  the  important  frequencies  while  rejecting  high-frequency  noise  in  the  magnitude 
spectrum  of  the  Tc  data.  Subsequently,  we  set  the  AR  model  order  m  =  5,  the  number  of  lags  in  the  data 
beyond  which  the  partial  autocorrelation  function  was  essentially  zero. 

We  used  the  measured  Tc  data  from  a  subject  to  estimate  b  using  the  standard  forward-backward  least 
squares  method  (see  [8],  chapter  8)  implemented  in  MATLAB  version  7.14  (function  ar).  We  also  estimated 
the  prediction  intervals  (Pis)  that  provide  information  about  the  uncertainty  of  the  predicted  values.  To 
estimate  the  Pis,  we  used  a  statistical  bootstrap  method,  where  a  population  of  models  is  built  based  on 
blocks  of  data  that  are  randomly  drawn  from  the  original  time  series  to  form  an  empirical  distribution  of 
models  (i.e.,  a  distribution  of  the  model  coefficients)  [9].  Following  this  procedure,  we  estimated  the 
covariance  matrix  Z  of  the  AR  model  from  a  distribution  of  models  for  an  M  time-step-ahead  predictor,  and 
used  the  following  equation  to  estimate  the  Pis  [10]: 

PI  =Za/2JyTZy  +  cj2  ,  (1) 

where  Za/2  denotes  the  prediction  factor  associated  with  an  a%  type  I  error,  y  represents  a  vector  of  data 
samples  y  -  \TCt  TCtA  .  .  .  7cf_m+1]r,  and  a"  denotes  the  variance  of  the  measurement  noise.  For  detailed 
information  about  the  AR  model  and  computation  of  the  Pis,  please  refer  to  Ref.  6. 

2.1.3  Sequential  Probability  Ratio  Test  Algorithm 

The  sequential  probability  ratio  test  (SPRT)  [11]  is  a  Bayesian  approach  that  considers  increasing  evidence 
from  a  sequence  of  observations  to  decide  whether  Tc  will  rise  beyond  a  given  temperature  threshold  and,  if 
so,  provides  an  alert  [12].  Briefly,  given  a  sequence  of  core  temperature  samples  Xu  X2,  ...  not  necessarily 
independent,  so  that  X  ~  5V  (ju^  ox)  is  a  normal  Gaussian  process  with  unknown  mean  px  and  a  given 
variance  cry2,  the  SPRT  tests  the  null  hypothesis  (H0)  that  /jx=  pQ  against  an  alternative  hypothesis  (Hi)  that 
pY  =  nv  where  pa  and  /.q  denote  the  mean  temperature  values  below  and  above  the  temperature  threshold, 
respectively,  with  pQ  <  /. i  If  p0  and  px  are  the  probability  density  functions  governing  H0  and  Hh 
respectively,  then  the  observed  likelihood  ratio  at  decision  time  t  can  be  represented  as  l,  =4A  p,(X,  k)  , 
where  K  is  the  length  of  the  sequence  of  samples  being  considered.  P0{X,  ,.) 
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In  order  to  apply  the  SPRT  algorithm  to  our  problem,  we  combined  three  predicted  values,  namely,  the  AR- 
model  prediction  (TCt_k),  the  upper  PI  (TCt_k  +  Pi,  k ),  and  the  lower  PI  ( TCt_k  -  Pi,  ^ ),  using  weights  0 
and  (/)  to  form: 

x,-k  =fCt_k--pit_ka-et-0),  (2) 

where  k  =  (0,  1,  K-  1)  denotes  a  time  index  and  the  weights  0 and  (j)  are  constrained  to  be  between  0  and 
1.  In  Eq.  (2),  when  (j)  -  0,  X,_k  equals  the  lower  PI;  when  <j)  —  I  and  0  -  0,  X,_k  equals  the  AR-model 
prediction;  and  when  0  =  <j)=  1,  X,.k  equals  the  upper  PI.  Thus,  X,  lies  between  the  lower  and  upper  Pis  for  all 
values  of  0  and  <j).  Note  that  a  decision  at  time  t,  for  a  prediction  window  of  P  min,  is  actually  made  at  time  t- 
P.  Then,  following  Wald’s  SPRT  methodology  [1 1],  we: 

accepted  H0  (below  temperature  threshold),  if  log(/,)  <  log(B);  or 

accepted  PL  (above  temperature  threshold),  if  log(/,)  <  log(A);  or 

made  no  decision  and  proceeded  to  time  t+\,  if  log(R)  <  log(/,)  <  log(A), 

where  A  and  B  are  constants  that  control  the  false-positive  rate  and  false-negative  rate,  respectively,  with  0  < 
B  <  A  <  co.  The  SPRT  algorithm  required  the  estimation  of  seven  parameters,  the  two  weights  0and  <j)  in  Eq. 

(2)  and  the  five  parameters  p(|,  crr  A,  and  B  from  a  subject’s  core  temperature  data.  For  detailed 
information  about  the  development  of  the  SPRT  algorithm,  please  refer  to  Ref.  6. 

2.1.4  Evaluating  the  Performance  of  the  Alert  Algorithm 

To  evaluate  the  alert  algorithm,  we  defined  an  “event”  as  an  episode  where  the  Tc  measurement  rises  and 
remains  above  a  specified  temperature  threshold  for  >  15  min.  The  event  ends  when  the  measured 
temperature  decreases  below  the  threshold  and  remains  below  the  threshold  for  >  15  min.  Thus,  an  event  was 
mapped  into  a  1  (true  response)  when  the  measured  temperature  was  above  the  threshold  and  0  otherwise. 
Similarly,  a  model-predicted  event  was  defined  as  an  episode  where  the  alert  algorithm’s  output  (the  model- 
predicted  response)  was  1 . 

We  evaluated  the  algorithm’s  performance  using  four  measures: 

(1)  Sensitivity:  the  fraction  of  time  points  during  which  both  the  true  and  model -predicted  responses  were  1. 
The  fraction  of  time  points  was  based  on  the  entire  time  series  of  Tc  measurements. 

(2)  Specificity:  the  fraction  of  time  points  during  which  both  the  true  and  model-predicted  responses  were  0. 
The  fraction  of  time  points  was  based  on  the  entire  time  series  of  Tc  measurements. 

(3)  Effective  prediction  horizon:  for  each  true  event,  we  computed  the  effective  prediction  horizon  by  adding 
the  prediction  window  P  to  the  time  difference  At  between  the  onset  of  the  true  event  and  the  onset  of  the 
model -predicted  event.  The  reported  effective  prediction  horizon  was  averaged  over  the  number  of 
events. 

(4)  Number  of  decision  switches:  the  cumulative  number  of  times  the  model-predicted  output  transitioned 
from  one  state  to  another  (0  to  1  or  1  to  0)  that  was  incongruent  with  the  true  state  (0  or  1)  of  the 
measured  temperature. 

We  computed  specificity  and  the  number  of  decision  switches  regardless  of  whether  or  not  an  event  had 
occurred,  however,  we  only  computed  sensitivity  and  effective  prediction  horizon  when  a  true  event 
occurred.  Note  that  while  specificity  measured  the  time  period  for  which  the  algorithm  incorrectly  predicted 
the  occurrence  of  an  event,  the  number  of  decision  switches  provided  the  number  of  times  the  algorithm 
predictions  incorrectly  switched  from  one  state  to  another. 
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We  evaluated  the  performance  of  the  alert  algorithm  on  the  six  subjects  whose  Tc  exceeded  38.5°C  using  a 
cross-validation  approach.  For  this  purpose,  we  estimated  the  AR  model  coefficients  b  and  the  SPRT 
parameters  on  one  subject’s  Tc  data,  applied  the  algorithm  on  the  other  five  subjects,  and  computed  the  four 
measures  of  performance  for  each  of  the  subjects.  We  repeated  this  procedure  for  each  of  the  six  subjects. 

2.2  Results:  Tc  Prediction  and  Alert  Algorithms 

Figure  2  shows  a  comparison  of  the  performance  of  the  SPRT -based  alert  algorithm  to  a  38.5°C  threshold 
applied  to  one  subject  (subject  #1).  In  Figure  2,  the  top  panel  shows  the  measured  data  (black  dash-dot  line), 
the  20-min-ahead  AR  model  predictions  (red  dashed  line),  and  the  corresponding  Pis  (shaded  region).  The 
middle  panel  depicts  the  true  events  (black  solid  line)  and  the  events  predicted  by  Prediction+PI  (red  dashed 
line),  while  the  bottom  panel  shows  the  same  plots  for  the  events  predicted  by  the  SPRT-based  alert 
algorithm.  Note  that  in  this  example,  the  parameters  of  the  AR  model  and  the  SPRT  algorithm  were 
estimated  using  the  Tc  data  from  subject  #6.  For  this  example,  Prediction+PI  yielded  100%  sensitivity,  95% 
specificity,  four  decision  switches  due  to  noisy  predictions,  and  an  effective  prediction  horizon  of  33  min 
(the  predicted  event  anticipated  the  true  event  by  13  min  beyond  the  20-min  prediction  window).  SPRT 
yielded  a  98%  sensitivity,  a  96%  specificity,  no  decision  switches,  and  an  effective  prediction  horizon  of  28 
min.  Thus,  SPRT  yielded  a  slightly  lower  sensitivity  and  a  reduced  effective  prediction  horizon  but, 
importantly,  produced  no  decision  switches  compared  with  Prediction+PI. 
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Figure  2:  The  top  panel  shows  the  measured  core  temperature  data  from  one  subject,  20-min-ahead  AR 
model  predictions,  and  the  corresponding  Pis.  The  middle  and  bottom  panels  show  the  depiction  of  a  true 
event  (solid  black  line)  and  algorithm  decisions  (dashed  red  line)  for  Prediction+PI  and  sequential  probability 
ratio  test  (SPRT)  algorithms,  respectively,  for  a  temperature  threshold  of  38.5°C. 
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Table  1 :  Performance  of  the  sequential  probability  ratio  test  (SPRT)  algorithm  applied  to  20-min-ahead 
auto-regressive  model  core  body  temperature  predictions  at  a  temperature  threshold  of  38.5°C 


Algorithm 

Sensitivity 

(%) 

Specificity 

(%) 

Effective 
prediction 
horizon  (min) 

Number  of 
decision 
switches 

Prediction+PI 

100.0  (0.0) 

70.6  (4.2) 

33  (3) 

14 

SPRT 

87.5  (7.6) 

92.2  (8.3) 

18(7) 

4 

Values  are  means  (standard  deviations)  from  a  fivefold  cross-validation  study  of  six 
subjects  whose  core  body  temperatures  exceeded  38.5°C  in  1 1  distinct  events. 


Table  1  shows  a  comparison  of  the  average  performance  of  the  SPRT-based  alert  algorithm  against  the 
Prediction+PI  algorithm  across  the  six  subjects  using  the  cross-validation  approach  described  in  Section 
2.1.4.  The  results  suggest  that  SPRT  increased  the  specificity  and  reduced  the  number  of  decision  switches  at 
the  cost  of  a  reduced  sensitivity  and  effective  prediction  horizon. 

2.3  Discussion:  Tc  Prediction  and  Alert  Algorithms 

While  Prediction+PI  yielded  the  largest  sensitivity  and  the  longest  effective  prediction  horizon,  it  also 
yielded  the  largest  number  of  decision  switches.  In  contrast,  delaying  alert  decisions  through  accumulation  of 
evidence  (SPRT)  yielded  a  small  number  of  decision  switches  at  the  cost  of  reduced  sensitivity  and  effective 
prediction  horizon.  For  practical  applications,  the  18-min  effective  prediction  horizon  yielded  by  SPRT 
provides  sufficient  time  to  enable  preventive  interventions  to  avoid  the  risk  of  heat  injuries.  Importantly,  the 
SPRT  algorithm  yielded  an  increased  specificity  and  a  significant  reduction  in  the  number  of  decision 
switches  and,  hence,  is  the  algorithm  recommended  for  the  development  of  a  real-time  early  warning  system. 

3.0  INDIVIDUALIZED  REAL-TIME  Tc  ESTIMATION  USING  NON-INVASIVE 
MEASUREMENTS 

As  mentioned  above,  a  practical  limitation  in  developing  an  early  warning  system  is  the  use  of  an  invasive 
sensor  to  measure  Tc.  To  obviate  the  necessity  of  an  ingestible  pill,  we  developed  a  mathematical  model  that 
uses  non-invasive  measurements  of  other  physiological  variables  and  two  environmental  variables  to  provide 
individualized  real-time  Tc  estimates.  We  detail  the  model  below. 

3.1  Methods:  Tc  Estimation  Using  Non-Invasive  Measurements 

3.1.1  Individualized  Tc  Estimation  Model 

The  proposed  individualized  model  uses  non-invasive  measurements  of  a  subject’s  Ac,  HR,  and  Ts  as  well  as 
two  environmental  variables,  TA  and  RH,  to  estimate  the  subject’s  T(  in  real  time  (Figure  3).  The 
individualized  model  includes  two  elements:  1)  a  mathematical  model  and  2)  a  Kalman  filter  [13].  First,  the 
mathematical  model  uses  the  measured  Ac  and  environmental  variables  TA  and  Ru  to  estimate  the  state 
variables  HR,  Ts,  and  Tc.  Then,  the  Kalman  filter  considers  the  error  between  the  measured  and  model- 
estimated  Hr  and  Ts  to  correct  the  state  variables  and  update  the  model  parameters,  resulting  in  an  improved 
estimation  and  individualization  of  Tc.  We  used  data  from  the  six  subjects  described  in  Section  2.1.1  for 
model  development  and  performance  evaluation. 

3.1.2  The  Mathematical  Model 

The  mathematical  model  consists  of  two  sub-models,  a  phenomenological  component  that  relates  Ac  to  HR 
and  a  first-principles,  macroscopic  energy  balance  component  [14]  that  relates  the  metabolic  activity 
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Figure  3:  The  individualized  model  for  core  body  temperature  (Tc).  The  inputs  to  the  model  are  the  measured 
heart  rate,  skin  temperature,  and  activity  profiles  from  a  subject  and  two  environmental  variables  (ambient 
temperature  and  relative  humidity).  The  activity  profile  and  environmental  factors  drive  the 
phenomenological  and  first-principles  models.  The  Kalman  filter  algorithm  then  considers  the  error  between 
the  measured  and  model-estimated  heart  rate  and  skin  temperature  to  update  the  model  parameters  and 

provide  individualized  real-time  Tc  estimates. 

(represented  by  HR )  to  Ts  and  Tc.  We  computed  Ac  from  triaxial  accelerometer  data  provided  through  the 
PSM  system  by  combining  the  three  orthogonal  axes  ( X ,  Y,  and  Z)  data  to  obtain  physical  activity  intensity. 
Then,  we  transformed  the  intensity  to  MET  values,  which  is  the  ratio  of  oxygen  consumed  during  physical 
activity  to  that  at  rest.  Next,  following  the  work  by  Sasaki  et  al.,  we  quantized  the  MET  values  into  five 
activity  levels:  0  for  rest  (MET  =  1),  1  for  light  activity  (MET  ~  1-3),  2  for  moderate  activity  (MET  ~  3-6),  3 
for  hard  activity  (MET  ~  6-9),  and  4  for  very  hard  activity  (MET  >  9)  [15].  Accordingly,  Ac  was  indexed  to 
represent  values  between  0  and  4.  The  phenomenological  component  was  represented  by  a  linear  first-order 
ordinary  differential  equation  (ODE)  driven  by  Ac  that  estimated  HR. 

The  first-principles  component  consists  of  two  ODEs.  The  first  ODE  represents  the  change  in  Tc  over  time 
due  to  heat  gained  in  the  body  from  metabolic  activity  (input  HR )  and  heat  lost  to  the  skin  due  to  blood  flow. 
The  second  ODE  represents  the  change  in  Ts  over  time  due  to  heat  gained  in  the  skin  from  the  body  and  heat 
lost  to  the  environment  via  the  difference  between  Ts  and  TA  and  sweat  evaporation  (which  depends  on  TA, 
Rh,  and  individual  factors,  such  as  clothing,  etc.).  The  mathematical  model  consists  of  three  ODEs  with 
seven  parameters. 

The  first  step  is  to  determine  a  set  of  initial  parameters  and  their  corresponding  SDs.  For  this  purpose,  we 
estimated  the  initial  model  parameters  using  data  from  all  subjects  (except  for  the  test  subject  for  whom  we 
wished  to  estimate  Tc )  by  fitting  the  model  to  these  data  so  as  to  minimize  the  overall  mean  squared  error 
between  the  model-estimated  physiological  variables  HR,  Ts,  and  Tc  and  their  corresponding  measurements. 
Subsequently,  we  estimated  the  parameters’  SDs.  We  repeated  this  procedure  for  each  of  the  six  subjects 
used  in  our  study,  thus  obtaining  six  slightly  different  “group-average”  models.  We  confirmed  that  all 
estimated  parameters  were  within  the  expected  range  obtained  from  experimental  studies  and  in  agreement 
with  previously  published  more  detailed  thermoregulatory  models  [14,  16,  17]. 

3.1.3  The  Kalman  Filter  Algorithm 

To  implement  the  Kalman  filter  algorithm,  we  converted  the  three  ODEs  of  the  mathematical  model  to  a 
discrete  time-varying  linear  model  with  the  model  states  consisting  of  HR,  Ts,  Tc,  and  the  model  parameters. 
We  initialized  the  Kalman  filter  with  the  group-average  model  parameters,  their  corresponding  SDs,  the 
noise  variances  of  the  HR  and  Ts  measurements  obtained  from  the  device  specifications,  and  a  measure  of  the 
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systemic  uncertainty  between  the  mathematical  model  estimates  and  a  subject’s  data.  We  estimated  the 
systemic  uncertainty  by  using  the  first  5  min  of  the  subject’s  data. 

After  initialization,  the  Kalman  filter  algorithm  proceeded  in  the  following  manner.  At  each  15-s  discrete 
time  interval,  the  algorithm  used  the  provided  TA  and  RH,  a  5-min  moving  average  of  the  measured  Ac  at  the 
current  time  point,  and  the  model  parameters  estimated  up  to  the  previous  time  point  to  update  the  model 
states.  Next,  the  algorithm  used  the  error  between  the  measured  HR  and  Ts  and  their  corresponding  model 
estimates  to  update  the  model  states  and  parameters  at  the  current  time  point.  The  algorithm  repeated  this 
procedure  for  each  time-step  until  the  end  of  the  time-series  data. 

3.1.4  Evaluation  of  the  Individualized  Model 

We  evaluated  two  variants  of  the  individualized  model:  1 )  one  variant  that  used  non-invasive  measurements 
of  Ac,  Hr,  and  Ts  along  with  hourly  measurements  of  the  environmental  variables  TA  and  RH  (original  model), 
and  2)  another  variant  that  only  used  the  measured  Ac  and  HR  and  month-average  values  of  the 
environmental  variables  (TA  =  27.2°C  and  RH  =  77%;  reduced  model).  The  reduced  model  allowed  us  to 
evaluate  model  performance  in  cases  where  Ts  measurements  and  hour-by-hour  weather  information  are  not 
available. 

As  discussed  in  above,  we  initialized  the  individualized  model  using  a  group-average  parameter  set  obtained 
while  excluding  the  test  subject  for  whom  we  wished  to  estimate  Tc.  This  allowed  us  to  cross  validate  the 
model  by  comparing  the  measured  and  estimated  Tc.  We  repeated  the  same  cross-validation  procedure  six 
times  to  assess  model  performance  for  each  of  the  six  subjects.  We  used  two  metrics  to  assess  the  model 
estimations:  the  root  mean  squared  error  (RMSE)  and  the  time  delay  between  the  measured  and  the  estimated 
Tc.  We  computed  the  RMSE  as  the  square  root  of  the  mean  squared  differences  between  the  model-estimated 
Tc  and  the  measured  Tc  across  the  time  duration  of  the  data.  We  computed  the  delay  as  the  time  difference 
between  the  maximum  values  of  the  auto-correlation  of  the  measured  Tc  and  the  cross-correlation  between 
the  measured  and  the  estimated  Tc. 

3.2  Results:  Tc  Estimation  Using  Non-invasive  Measurements 

Figure  4  shows  the  performances  of  the  two  variants  of  the  individualized  model  on  one  subject  (subject  #1). 
The  top-left  panel  shows  the  measured  Ac  (gray  dash-dot  line)  and  its  5-min  moving  average  (solid  black 
line).  The  Ac  drove  the  individualized  model  to  track  the  measured  HR  and  Ts  (middle  and  bottom  panels, 
respectively)  via  the  Kalman  filter,  and  provided  real-time  Tc  estimates  (right  panel).  The  original  and 
reduced  models  estimated  Tc  with  RMSEs  of  0.2 1°C  and  0.24°C  and  delays  of  0  min  and  1  min, 
respectively. 

Table  2  shows  the  average  RMSEs  and  estimation  delays  for  the  original  and  reduced  models  across  the  six 
subjects.  The  original  and  reduced  models  yielded  average  RMSEs  of  0.3 1°C  (SD  0.09)  and  0.35°C  (SD 
0.10)  and  average  estimation  delays  of  3  min  (SD  4)  and  5  min  (SD  4),  respectively.  Thus,  the  reduced  model 
yielded  a  13%  higher  RMSE  and  a  2-min  additional  delay  compared  with  the  original  model. 

3.3  Discussion:  Tc  Estimation  Using  Non-invasive  Measurements 

A  novel  feature  of  the  individualized  model  is  that  it  adapts  the  model  parameters  to  a  subject  and,  as  such, 
explicitly  accounts  for  subject-specific  variations  in  thermoregulatory  responses,  acclimation,  and  responses 
to  exogenous  factors,  such  as  clothing  and  environmental  conditions.  For  example,  if  RH  increases  or  a 
subject’s  clothing  impedes  sweat  evaporation,  the  model  accounts  for  these  changes  by  reducing  the  rate  of 
heat  loss  due  to  sweat  evaporation,  leading  to  an  increase  in  Tc.  This  is  possible  because  the  Kalman  filter 
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Figure  4:  Performance  of  two  variants  of  the  individualized  model  on  one  subject  (subject  #1).  The  original 
model  used  the  measured  activity,  heart  rate  [in  beats/min  (bpm)],  and  skin  temperature  along  with  hourly 
measurements  of  two  environmental  variables:  ambient  temperature  and  relative  humidity.  The  reduced 
model  used  month-average  values  of  the  environmental  variables  and  did  not  use  the  measured  skin 
temperature.  The  top-left  panel  shows  the  measured  activity  (gray  dash-dot  line)  and  a  5-min  moving  average 
(solid  line).  The  middle-  and  bottom-left  panels  show  the  measured  heart  rate  and  skin  temperature, 
respectively,  and  the  corresponding  model  estimates.  The  right  panel  shows  the  measured  core  temperature 

and  the  corresponding  model  estimates. 


uses  the  non-invasive  measurements  of  HR  and  Ts  to  update  the  model  parameters  in  real  time  and  adapt  them 
to  subject-specific  responses. 


The  two  variants  of  the  individualized  model  estimated  Tc  with  errors  of  0.3 1°C  and  0.35°C  and  small  delays 
of  3-  to  5-min  delay  (Table  2).  These  errors  are  within  the  physiological  variations  of  Tc  in  an  individual 
(0.40°C)  [18].  The  estimation  delays  occurred  mainly  because  of  the  inherent  time  lags  associated  with  the 
low-pass  filters  used  to  smooth  the  HR  and  Ts  measurements.  While  an  “optimal”  algorithm  should  provide 
Tc  estimates  with  no  delay,  we  believe  that  such  a  requirement  may  not  be  practical  for  real-time  applications 

Table  2:  Performances  of  the  two  variants  of  the  individualized  models  across  the  six 
subjects.  The  original  model  estimated  Tc  using  non-invasive  measurements  of  Ac,  Hr,  and  Ts 
and  hourly  measurements  of  TA  and  Rh ■  The  reduced  model  estimated  Tc  using  non-invasive 
measurements  of  Ac  and  Hr  and  month-average  values  of  TA  and  Rh 


Individualized 

models 

RMSE  (°C) 

Delay  (min) 

Original  model 

0.31  (0.09) 

3(4) 

Reduced  model 

0.35  (0.10) 

5(4) 

Values  are  means  (standard  deviations).  Ac,  activity  level;  HR, 
heart  rate;  RH,  relative  humidity;  TA,  ambient  temperature;  Tc,  core 
body  temperature;  Ts,  skin  temperature. 
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when  the  time-series  signal  is  noisy  and  needs  to  be  further  smoothed  before  use  [5].  Moreover,  a  3-  to  5-min 
delay  represents  less  than  15%  of  the  time  taken  by  Tc  to  rise  by  1°C  during  physical  activity  [16].  Thus,  the 
Tc  estimates  provided  by  the  individualized  model  can  be  used  as  a  surrogate  for  the  ingestible  pill  and  drive 
the  alert  algorithm  detailed  in  Section  2.0.  This  could,  for  the  first  time,  allow  us  to  develop  a 
hardware/software  system  for  real-time  warning  of  an  impending  heat  injury. 

Very  few  of  the  existing  non-invasive  performance-  and  fitness-tracking  devices  provide  Ts  measurements, 
e.g.,  Basis  Peak™  (Intel  Inc.,  San  Francisco,  CA).  Hence,  the  observation  that  the  reduced  model  that  only 
requires  HR  and  Ac,  which  are  readily  available  in  many  devices,  was  only  13%  worse  than  the  original 
model  has  an  important  ramification  in  the  development  of  a  reliable  real-time  warning  system. 

4.0  CONCLUSIONS 

In  conclusion,  we  developed  an  alert  algorithm  to  provide  reliable  ahead-of-time  warning  of  an  impending 
rise  in  Tc.  We  showed  that  20-min-ahead  AR  model  predictions  combined  with  the  SPRT  algorithm  provides 
a  sufficiently  long  effective  prediction  horizon  (18  min),  with  a  high  sensitivity  and  specificity  (>  87%)  and  a 
small  number  of  decision  switches  (four)  to  enable  preventive  interventions  to  minimize  the  risk  of  heat 
injuries.  We  then  developed  an  individualized  model  that  uses  non-invasive  measurements  of  physiological 
signals  and  environmental  variables  to  adapt  the  model  parameters  and  provide  subject-specific, 
individualized  Tc  estimates.  We  found  that  the  individualized  model  provided  accurate  Tc  estimates  even 
when  we  used  month-averaged  values  of  the  environmental  variables  and  neglected  skin  temperature 
measurements.  These  findings  suggest  that  the  model-estimated  Tc  can  replace  the  invasive  core  temperature 
measuring  pill  and  thus  enable  the  development  of  an  early  warning  system  that  can  be  deployed  in 
ambulatory  settings.  Currently,  we  are  in  the  process  of  integrating  this  model  with  the  alert  algorithm  in  a 
smartphone  application  to  provide  reliable  early  warning  of  an  impending  rise  in  core  temperature,  which 
would  help  prevent  heat  injuries. 
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